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Abstract 

' This paper is concerned with the analysis of convergent sequential and parallel over- 

, lapping domain decomposition methods for the minimization of functionals formed by a 

discrepancy term with respect to data and a total variation constraint. To our knowl- 
edge, this is the first successful attempt of addressing such strategy for the nonlinear, 
nonadditive, and nonsmooth problem of total variation minimization. We provide several 
numerical experiments, showing the successful application of the algorithm for the restora- 
tion of ID signals and 2D images in interpolation/inpainting problems respectively, and 
in a compressed sensing problem, for recovering piecewisc constant medical-type images 
^ ' from partial Fourier ensembles. 
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^ ! 1 Introduction 



In concrete apphcations, e.g., for image processing, one might be interested to recover at best 
a digital image provided only partial linear or nonlinear measurements, possibly corrupted 
by noise. Given the observation that natural and man-made images can be characterized by 
a relatively small number of edges and extensive relatively uniform parts, one may want to 
help the reconstruction by imposing that the interesting solution is the one which matches 
the given data and has also a few discontinuities localized on sets of lower dimension. 

In the context of compressed sensing O [TJ [HI [21], it has been clarified that the min- 
imization of £i-norms occupies a fundamental role for the promotion of sparse solutions. 
This understanding furnishes an important interpretation of total variation minimization, 
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i.e., the minimization of the L^-norm of derivatives [3l], as a regularization technique for 
image restoration. The problem can be modelled as follows; let Q, C M'^, for d = 1,2 be a 
bounded open set with Lipschitz boundary, and H = L^(0). For u € Lj^^{Q) 

sup ^ u div (p dx : (p G [C^(r2)] , H'/'lloo < 1 > 

is the variation of u. Further, u G BV{Q), the space of bounded variation functions [HIM]; if 
and only if V{u,Q) < oo. In this case, we denote \D{u)\{^) = V{u,Q). If w € W^'^{^1) 
(the Sobolev space of L^-functions with L^-distributional derivatives), then |L'(n)|(r2) = 
|Vn| dx. We consider as in [12\ [38] the minimization in BV{Q) of the functional 

J{u) :=\\Tu- g\\l + 2a\D{u)\{n), (1) 

where T : L^(r2) — > L^(0) is a bounded linear operator, g € L^(il) is a datum, and a > 
is a fixed regularization parameter [23]. Several numerical strategies to perform efficiently 
total variation minimization have been proposed in the literature. Without claiming of being 
exhaustive, we list a few of the relevant methods, ordered by their chronological appearance: 

(i) the linearization approach of Vogel et al. [20] and of Chambolle and Lions [12] by 
iteratively re-weighted least squares, see also [18] for generalizations and refinements in the 
context of compressed sensing; 

(ii) the primal-dual approach of Chan et al. [13] : 

(iii) variational approximation via locally quadratic functionals as in the work of Vese et 
al. [21 [38]; 

(iv) iterative thresholding algorithms based on projections onto convex sets as in the work 
of Chambolle [10] as well as in the work of Combettes and Wajs [H] and Daubechies et al. 

CS]; 

(v) iterative minimization of the Bregman distance as in the work of Osher et al. [33] 
(also notice the very recent Bregman split approach [27]): 

(vi) graph cuts [11[ [T6] for the minimization of ([1]) with T = Id and an anisotropic total 
variation; 

(vii) the approach proposed by Nesterov [31] and its modifications by Weiss et al. [39| . 
These approaches differ significantly, and they provide a convincing view of the interest 

this problem has been able to generate and of his applicative impact. However, because of 
their iterative-sequential formulation, none of the mentioned methods is able to address in 
real-time, or at least in an acceptable computational time, extremely large problems, such 
as 4D imaging (spatial plus temporal dimensions) from functional magnetic-resonance in 
nuclear medical imaging, astronomical imaging or global terrestrial seismic tomography. For 
such large scale simulations we need to address methods which allow us to reduce the problem 
to a finite sequence of sub-problems of a more manageable size, perhaps computable by one 
of the methods listed above. With this aim we introduced subspace correction and domain 
decomposition methods both for ^i-norm and total variation minimizations [25 [ I26 [ [35]. We 
address the interested reader to the broad literature included in [26] for an introduction to 
domain decompositions methods both for PDEs and convex minimization. 

1.1 Difficulty of the problem 

Due to the nonsmoothness and nonadditivity of the total variation with respect to a nonover- 
lapping domain decomposition (note that the total variation of a function on the whole domain 
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equals the sum of the total variations on the subdomains plus the size of the jumps at the 
interfaces |26l formula (3.4)]), one encounters additional difficulties in showing convergence 
of such decomposition strategies to global minimizers. In particular, we stress very clearly 
that well-known approaches as in [9l [HI |36l EZ] are not directly applicable to this problem, 
because either they do address additive problems or smooth convex minimizations, which 
is not the case of total variation minimization. Moreover the interesting solutions may be 
discontinuous, e.g., along curves in 2D. These discontinuities may cross the interfaces of the 
domain decomposition patches. Hence, the crucial difficulty is the correct numerical treat- 
ment of interfaces, with the preservation of crossing discontinuities and the correct matching 
where the solution is continuous instead, see [26l Section 7.1.1]. 

The work [26] was particularly addressed to nonoverlapping domain decompositions f^i U 
^^2 C C U and Pi ^2 = 0. Associated to the decomposition define = {n € 
L^(r2) : supp(n) C f^j}, for i = 1,2; note that L'^{U) = Vi © V2. With this splitting we wanted 
to minimize by suitable instances of the following alternating algorithm: Pick an initial 
Vi (BV2 3 + u^2^ := u^^\ for example u'^^^ = 0, and iterate 

f (n+l) • , (n.)N 

u\ PS argmm^,^gVi J{Vi + ^ ') 

(n+l) _ . ^/ (n+l) . X 

< ~ argmmt,2gv'2 ■^(^1 + ^2) 

^ ^(n+i) :=^(-+i)+4«+i). 

In [26j we proposed an implementation of this algorithm which is guaranteed to converge 
and to decrease the objective energy J monotonically. We could prove its convergence to 
minimizers of J only under technical conditions on the interfaces of the subdomains. However, 
in our numerical experiments, the algorithm seems always converging robustly to the expected 
minimizer. This discrepancy between theoretical analysis and numerical evidences motivated 
our investigation on overlapping domain decompositions. The hope was that the redundancy 
given by overlapping patches and the avoidance of boundary interfaces could allow for a 
technically easier theoretical analysis. 

1.2 Our approach, results, and technical issues 

In this paper we show how to adapt our previous algorithm [26] to the case of an overlap- 
ping domain decomposition. The setting of an overlapping domain decomposition eventually 
provides us with a framework in which we successfully prove its convergence to minimizers of 
J^, both in its sequential and parallel forms. Let us stress that to our knowledge this is the 
first method which addresses a domain decomposition strategy for total variation minimiza- 
tion with a formal theoretical justification of convergence. It is important to mention that 
there are other very recent attempts of addressing domain decomposition methods for total 
variation minimization with successful numerical results |30j . 

Our analysis is performed for a discrete approximation of the continuous functional ([T|), 
for ease again denoted J" in ([3]). Essentially we approximate functions u by their sampling 
on a regular grid and their gradient Du by finite differences Vu. It is well-known that such 
a discrete approximation F-converges to the continuous functional (see [4]). In particular, 
discrete minimizers of ([3]), interpolated by piecewise linear functions, converge in weak-*- 
topology of BV to minimizers of the functional ^ in the continuous setting. Of course, 
when dealing with numerical solutions, only the discrete approach matters together with its 
approximation properties to the continuous problem. However, the need of working in the 
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discrete setting is not only practical, it is also topological. In fact bounded sets in BV are 
(only) weakly-*-compact, and this property is fundamental for showing that certain sequences 
have converging subsequences. Unfortunately, the weak-*-topology of BV is "too weak" for 
our purpose of proving convergence of the domain decomposition algorithm; for instance, 
the trace on boundary sets is not a continuous operator with respect to this topology. This 
difficulty can be avoided, for instance, by F-approximating the functional ([1]) by means of 
quadratic functional (as in [2l|T2l[38]) and working with the topology of ^^^'^(17), the Sobolev 
space of L^-functions with L^-distributional first derivatives. However, this strategy changes 
the singular nature of the problem which makes it both interesting and difficult. Hence, 
the discrete approach has the virtues of being practical for numerical implementations, of 
correctly approximating the continuous setting, and of retaining the major features which 
makes the problem interesting. Note further that in the discrete setting where topological 
issues are not a concern anymore, also the dimension d can be arbitrary, contrary to the 
continuous setting where the dimension d has to be linked to boundedness properties of the 
operator T, see [38l property H2, pag. 134]. For ease of presentation, and in order to avoid 
unnecessary technicalities, we limit our analysis to splitting the problem into two subdomains 
r^i and r22- This is by no means a restriction. The generalization to multiple domains comes 
quite natural in our specific setting, see also [26l Remark 5.3]. When dealing with discrete 
subdomains Oj, for technical reasons, we will require a certain splitting property for the total 
variation, i.e., 

\Vu\{n) = |Vn|nJ(l^i) + ci(u|(n2\ni)uri), \Vu\{n) = |Vu|n2 |(f^2) + C2(n|(f^^\f^2)ur2)> (2) 

where ci and C2 are suitable functions which depend only on the restrictions 'u|(f72\f^i)uri 
and u|(r2i\Q2)ur2 respectively, see ([9]) (symbols and notations are clarified once for all in the 
following section). Note that this formula is the discrete analogous of [26\ formula (3.4)] 
in the continuous setting. The simplest examples of discrete domains with such a property 
are discrete d-dimensional rectangles (d-orthotopes). Hence, for ease of presentation, we will 
assume to work with d-orthotope domains, also noting that such decompositions are already 
sufficient for any practical use in image processing, and stressing that the results can be 
generalized also to subdomains with different shapes as long as ([2]) is satisfied. 

1.3 Organization of the work 

The paper is organized as follows. In Section [2] we collect the relevant notations and symbols 
for the paper. Section [3] introduces the problem and the overlapping domain decomposition 
algorithm which we want to analyze. In Section U] we address the approximate solution of the 
local problems defined on the subdomains Qi and we show how we interface them, by means 
of a suitable Lagrange multiplier. Section [5] and Section [6] are concerned with the convergence 
of the sequential and parallel forms of the algorithm. In particular, in Section [5] we provide a 
characterization of minimizers by a discrete representation of the subdifferential of J^. This 
characterization is used in the convergence proofs in order to check the reached optimality. 
The final Section [7] provides a collection of applications and numerical examples. 

2 Notations 

Let us fix the main notations. Since we are interested in a discrete setting we define the 
discrete d-orthotope fl = {x\ < . . . < x]^ } x . . . x {xf < . . . < xf^ } C M"^, d G N and the 
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considered function spaces are H = ]^^ix^2x...x7Vrf^ where A^j S N for i = 1, . . . , d. For u £ Ti 
we write u = n(xi)igi with 

d 

I:=ll{l,...,Nk} 

k=l 

and 

u{xi) =u{xl,...,xfj 
where € {1, . . . , N^} and (xi)igj € ft. Then we endow Ti with the norm 

1/2 / \ 1/2 

2 I 



We define the scalar product of u,v G H as 



and the scalar product of p, g € W'' as 

with (y, = J2j=i VjZj for every y = (yi, . . . , y^) e K'^ and z = (zi, . . . , z^) G M"^. We will 
consider also other norms, in particular 

i/p 



(Y.\u{xi)\P^ , l<p<oo, 



and 

llulloo = SUp|'u(Xi)|. 
iGX 

We denote the discrete gradient Vu by 

(Vn)(xi) = ((Vn)Hxi),...,(V^x)'^(xi)) 

with 

{VuYixi) = h^""^' ' ■ ■ ■ ' ^^^+1' • • • ' ^^<^) " ""^^^'i ' • • • ' 4 ' • • • ' 4) if < 

for all J = 1, . . . ,d and for all i = (ii, . . . , i^) € T. 
Let V3 : M ^ M, we define for u £ 



^{\u;\m =Y^i\u;ixO\) = Y^iHx)\), 



where |y| = y y^ + . . . + y^. In particular we define the total variation of u by setting (p{s) = s 
and (jj = Vu, i.e., 

\\/u\{n) := \Vu{xi)\ = Y |Vn(x)|. 
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For an operator T we denote T* its adjoint. Further we introduce the discrete divergence 
div : T-C^ —> H defined, in analogy with the continuous setting, by div = —V* (V* is the 
adjoint of the gradient V). The discrete divergence operator is exphcitly given by 



{divp){xi) 



'p'^{xl,...,xfj 
p'^{xl,...,xfj 



...+ < 



iil<ii<Ni 
if ii = 1 
if n = iVi 

ifl<id<Nd 
if id = 1 
if id = Nd, 



for every p = {p^,...,p'^) G T-C^ and for ah i = (ii, . . . , id) G I- (Note that if we considered 
discrete domains which are not discrete d-orthotopes, then the definitions of gradient and 
divergence operators should be adjusted accordingly.) With these notations, we define the 
closed convex set 

K := jdivp -.pen"^, \p{x)\^ < 1 for all x G f]} , 

where = max . . . , |p'^(x)|}, and denote Pk{u) = argmin^g^ ||n — ■u||2 the 

orthogonal projection onto K. We will often use the symbol 1 to indicate the constant vector 
with entry values 1 and 1_d to indicate the characteristic function of the domain C fi. 



3 The Overlapping Domain Decomposition Algorithm 

We are interested in the minimization of the functional 

J{u) := \\Tu -g\\l + 2a \V{u)\ (fi). 



(3) 



where T G jC{Tt) is a linear operator, 51 G is a datum, and a > is a fixed constant. In 
order to guarantee the existence of minimizers for ([3]) we assume that: 

(C) J is coercive in Tl, i.e., there exists a constant C > such that {J^ < C} := {u G : 
J^{u) < C} is bounded in H. 

It is well known that if 1 ^ ker(T) then condition (C) is satisfied, see \38\ Proposition 3.1]. 

Now, instead of minimizing ([3]) on the whole domain we decompose ft into two overlapping 
subdomains Qi and Q2 such that 17 = l^i U ^12, ^1 r\ 0,2 7^ 0, and ([2]) is fulfilled. For 
consistency of the definitions of gradient and divergence, we assume that also the subdomains 
rij are discrete (i-orthotopes as well as fi, stressing that this is by no means a restriction, but 
only for ease of presentation. Due to this domain decomposition Ti is split into two closed 
subspaces Vj = {u £ Ti : supp(n) C ilj}, for j = 1, 2. Note that H = Vi + V2 is not a direct 
sum of Vi and V2, but just a linear sum of subspaces. Thus any n G H has a nonunique 
representation 



u[x 



ui{x) X G r^i \ 

ui{x) + U2{x) X G $7i n S72 ) Ui£Vi, i 

^U2{x) X £ 0,2X^1 



1,2. 



(4) 
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We denote by Fi the interface between Qi and ^2 \ ^^i and by r2 the interface between Q2 
and \ (the interfaces are naturally defined in the discrete setting). 
We introduce the trace operator of the restriction to a boundary 

TV|r,:^i^Mr% i = l,2 

with Tr |r- Vi = Vi Ir,, the restriction of on Fj. Note that M^' is as usual the set of maps 
from Fj to M. The trace operator is clearly a linear and continuous operator. We additionally 
fix a bounded uniform partition of unity (BUPU) {xi,X2} C 7i such that 

(a) T\^|r, Xi = fori = 1,2, 

(b) Xi+X2 = 1, 

(c) suppxi C $7j for i = 1,2, 

(d) max{||xi||oo, IIX2II00} = k < 00. 
We would like to solve 

argmin„g^ J{u) 

by picking an initial Vi + V2 3 uf'^ + u^2^ := u^^^ S 7i, e.g., uf^^ = 0, i = 1, 2, and iterate 

(n+l) ■ rrl , ~{")\ 

u\ ' « argmm ^^^v^ J\v\^u\') 
Tr|rit'i=0 

u\ ^ argmm ^jgy^ ^(^1 + ""2) 

""l + ^2 



^(n+1) ._ ,,{«+!) _L 



~(n+l) 

u\ 

~(n+l) 



X2 



Note that we are minimizing over functions f j € for i = 1,2 which vanish on the interior 
boundaries, i.e., Tr jr^ Vi = 0. Moreover is the sum of the local minimizers and 

(n) 

U2 , which are not uniquely determined on the overlapping part. Therefore we introduced 
a suitable correction by xi and X2 in order to force the subminimizing sequences {u^^^)neN 
and (n2"'*)nGN to keep uniformly bounded. This issue will be explained in detail below, see 
Lemma 15.51 From the definition of x«; i = 1, 2, it is clear that 

Note that in general ti^"^ 7^ and tig"^ 7^ ^'2"^- -'^^ ® "~" (^'^^ approximation 

symbol) because in practice we never perform the exact minimization. In the following section 
we discuss how to realize the approximation to the individual subspace minimizations. 

4 Local Minimization by Lagrange Multipliers 

Let us consider, for example, the subspace minimization on 0,i 

argmin „^6Vi + n2) = argmin ^^^v^ \\Tvi - {g - Tu2)\\l + 2a\V{vi + U2 \ni)\ {^). 

Tr|r^?)i=0 Tr|rjfi=0 

(6) 
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First of all, observe that {u & TC : Tr {r^^ u = Tr |ri U2, J{u) < C} C {J^ < C}, hence the 
former set is also bounded by assumption (C) and the minimization problem ^ has solutions. 

It is useful to us to introduce an auxiliary functional J'j' of J^, called the surrogate func- 
tional of i7 (cf. [26]): Assume a,ui G Vi, € V2, and define 

Jt{ui+U2,a) := J{ui+U2) + ||ni - a\\l - \\T{ui - a)\\l (7) 

A straightforward computation shows that 

J{{ui + U2,a) = llm - (a + {T*{g - Tu2 - Ta)) + 2a \V{ui + U2)\ (0) + $(0,5,^2), 

where # is a function of a,g,U2 only. Note that now the variable ui is not anymore effected 
by the action of T. Consequently, we want to realize an approximate solution to ([6]) by using 
the following algorithm: For uf'^ = uf'^ G ^i, 

= argmin ^^gvi J{{ui U2,uf^), £>0. (8) 

Tr|riiii=0 

Additionally in ([8]) we can restrict the total variation on f^i only, since we have 

|V(ni +^2)1 (O) = |V('Ui + 7/2) Is^J +C2(n2|(n2\!^,)uri)- (9) 

where we used ([2]) and the assumption that ui vanishes on the interior boundary Fi. Hence 
dS]) is equivalent to 

argmin ^^^v^ Jf'(tti + n2, ^) = argmin ^^^v^ \\ui - zi\\l + 2a \V {ui -\- U2) \nj\ {^1) , 

Tr|ri«i=0 Tr|ri«i=0 

where zi = u^^ + (T* {g — Tu2 — Tu^p)) |qj . Similarly the same arguments work for the second 
subproblem. 

Before proving the convergence of this algorithm, we need to clarify first how to practically 
compute u^^"*"^^ for uf' given. To this end we need to introduce further notions and to recall 
some useful results. 

4.1 Generalized Lagrange multipliers for nonsmooth objective functions 

Let us begin this subsection with the notion of a subdifferential in finite dimensions. 

Definition 4.1. For a finite locally convex space V and for a convex function F : V ^ 
M U {—00, +00}, we define the subdifferential of F at x (zV, as the set valued function 



dF{x) :-- 



if F{x) = 00 

{x* S V : {x*, y — x) + F{x) < F{y) \/y E V} otherwise. 



It is obvious from this definition that E dF{x) if and only if x is a minimizer of F. Since we 
deal with several spaces, namely, Ti.,Vi, it will turn out to be useful to sometimes distinguish 
in which space the subdifferential is defined by imposing a subscript dyF for the subdifferential 
considered on the space V . 
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We consider the following problem 

argmin^gy{F(a;) : Gx = 6}, (10) 
where G : V ^ V is a linear operator on V. We have the following useful result. 

Theorem 4.2. JM, Theorem 2.1.4, p. 305] Let N = {G*\ : A G = Range(G*). Then, 
Xo S {x S y : G{x) = h} solves the constrained minimization problem (jlOp if and only if 

G dF{xQ) + N. 

4.2 Oblique thresholding 

We want to exploit Theorem 14.21 in order to produce an algorithmic solution to each iteration 
step dSI), which practically stems from the solution of a problem of this type 

argmin „jgvi \\ui - zi\\\ + 2a\'\I {ui + U2 

Tr|ri«i=0 

It is well-known how to solve this problem if U2 = in Cti and the trace condition is not 
imposed. For the general case we propose the following solution strategy. In what follows all 
the involved quantities are restricted to e.g., ui = ui |nn^^2 = U2 

Theorem 4.3 (Oblique thresholding). For U2 G V2 and for zi G Vi the following statements 
are equivalent: 

(i) u\ = aigmiu - zilH + 2a |V(ui + ^2)! (J^i); 

Tr|r^«i=0 

(a) there exists rj G Range(Tr |r^)* = {?? G Vi with supp(r/) = Fi} such that G n]^ — (zi — 
r^) + adv, |V(- + U2)|(0i)(nt); 

(Hi) there exists r/ G Vi with supp(?7) = Fi such that u\ = [I — PaK){zi + U2 — f]) — U2 ^Vi 
and Tr |ri ^1 = 0; 

(iv) there exists rj ^Vi with supp(?7) = Fi such that Tr Ir^ ?7 = Tr Ir^ zi + Tr Ir^ PaKij] — 
{zi + U2)) or equivalently rj = (Tr IrJ* Tr |ri {zi + PaKiv - {zi + U2))). 

We call the solution operation provided by this theorem an oblique thresholding, in analogy 
to the terminology in [T7j, because it performs a thresholding of the derivatives, i.e., it sets 
to zero most of the derivatives oi u = U1+U2 ^ zi on Oi, provided U2 which is a fixed vector 
in V2. 

Proof. Let us show the equivalence between (i) and (ii). The problem in (i) can be reformu- 
lated as 

ul = argmin„^gyjF(ui) := ||iii - zi\\l + 2a |V(ui + U2) \ (f^i),Tr |ri ui = 0}. (11) 

Recall that Tr jr^: Vi M^^ is a surjective map with closed range. This means that 
(Tr jpi)* is injective and that Range(Tr IrJ* = {rj £ Vi with supp(?7) = Fi} is closed. Using 
Theorem 14.21 the optimality of is equivalent to the existence of 77 G Range(Tr |ri)* such 
that 

0edv,F{ul) + 2ri. (12) 



OVERLAPPING DOMAIN DECOMPOSITION METHODS FOR TV-MINIMIZATION 10 



Due to the continuity of ||ui — in Vi, we have, by [221 Proposition 5.6], that 

dv,F{ul) = 2{ul - zi) + 2adv, |V(- + U2)\ (f^i)K). (13) 

Thus, the optimahty of is equivalent to 

Oeul-zi+rj + adv^ |V(- + U2)\ {^i){u[). (14) 

This concludes the equivalence of (i) and (ii). Let us show now that (iii) is equivalent to (ii). 
The condition in (iii) can be rewritten as 

^* = (/-P^^)(zi +U2 -r?), ^*=u\+U2. 

Since |V(-)| > is 1-homogeneous and lower-semicontinuous, by [261 Example 4.2.2], the 
latter is equivalent to 

G r - (^1 + ^2 - ??) + adv, |V(-)| i^iW). 

and equivalent to (ii). Note that in particular we have dy-^ |V(-)| (r2i)(^*) = dy-^ |V(- + U2)\ (r2i)(nj), 
which is easily shown by a direct computation from the definition of subdifferential. We prove 
now the equivalence between (iii) and (iv). We have 

u\ = {I - PaK){zi + U2 - ri) - U2 eVi, T] & Vi with supp(r/) = Ti, Tr |ri n]; = 

= Zi-7]- PaK{zi +U2- rf). 

By applying Tr Ir^ to both sides of the latter equality we get 

= Tr |ri zi - Tr |ri - Tr |ri PaK{zi + U2 - r]). 
By observing that — Tr PaxiO = Tr -Pa_ft'(— 0) we obtain the fixed point equation 

Tr |ri = Tr |ri zi + Tr |ri Paxiv - {zi + U2)). (15) 
Conversely, since all the considered quantities in 

(/ - PaK)izi +U2-r]) -U2 

are in Vi , the whole expression is an element in Vi and hence as defined in (iii) is an element 
in Vi and Tr |ri = 0. This shows the equivalence between (iii) and (iv) and therewith 
finishes the proof. □ 

We wonder now whether any of the conditions in Theorem 14.31 is indeed practically satis- 
fied. In particular, we want to show that i] G Vi as in (iii) or (iv) of the previous theorem is 
provided as the limit of the following iterative algorithm: 

r/W G Vi,supp7?W = Ti 7?("^+i) = (Tr |rJ*Tr \r, [zi + P^k {^^""^ - {zi + U2))) , m > 0. 

(16) 

Proposition 4.4. The following statements are equivalent: 

(i) there exists rj ^Vi such that t] = (Tr |ri)* Tr |ri {zi + Paxiv — (zi + ^^2))) (which is in 
turn the condition (iv) of Theorem \4.3\ ) 
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(ii) the iteration U6\} converges to any G Vi that satisfies U5\) . 

For the proof of this Proposition we need to recall some well-known notions and results. 

Definition 4.5. A nonexpansive map T : ?{ ^ Ti is strongly nonexpansive if for [un — Vn)n 
hounded and \\T{un) — T{vn)\\2 — W^n — Vn\\2 ^ we have 

Un- Vn- {T{Un) -T{Vn)) ^ 0, n OO. 

Proposition 4.6 (Corollaries 1.3, 1.4, and 1.5 [5])' LetT : Tl ^ 7i be a strongly nonexpansive 
map. Then hxT = {u € 7i : T{u) = u\ ^ % if and only if {T"'u)n converges to a fixed point 
uq € HxT for any choice of u € Ti.. 

Proof. (Proposition 14. 4( ) Projections onto convex sets are strongly nonexpansive [31 Corollary 
4.2.3]. Moreover, the composition of strongly nonexpansive maps is strongly nonexpansive 
[21 Lemma 2.1]. By an application of Proposition 14.61 we immediately have the result, since 
any map of the type T(,^) = Q{£,) + £,o is strongly nonexpansive whenever Q is (this is 
a simple observation from the definition of strongly nonexpansive maps). Indeed, we are 
looking for fixed points of r/ = (Tr |ri)* Tr Ir^ (^i + Paxiv — {zi + U2))) or, equivalently, of 
C = (Tr |rJ*TV \r, P„i^(e)-((TV |rJ*Tr |r, U2), where C = (TV |rJ*Tr \r, {v-{zi+U2)). □ 

^ V ' ^ V ' 

■=Q :=5o 

4.3 Convergence of the subspace minimization 

Prom the results of the previous section it follows that the iteration dSj) can be explicitly 
computed by 

nf = S^inf^ +T*{g- Tu2 - Tnf^) + U2 - v^'^) " ^2, (17) 
where Sa '■= I — PaK and G Vi is any solution of the fixed point equation 

7? = (TV IrJ* Tr |r, [{uf +T*{g- Tu^ - Tu'f)) - Paxiuf^ +T*{g- Tu2 - Tuf^ + U2 - v))) ■ 
The computation of rj^^^ can be implemented by the algorithm (jl6p . 

Proposition 4.7. Assume U2 S V2 and \\T\\ < 1. Then the iteration (|17|) converges to a 
solution ul G Vi 0/ dnj) for any initial choice of uf''' G Vi . 

The proof of this proposition is standard, see [T5l [T71 [26]. 

Let us conclude this section mentioning that all the results presented here hold symmetri- 
cally for the minimization on V2, and that the notations should be just adjusted accordingly. 

5 Convergence of the Sequential Alternating Subspace Mini- 
mization 

In this section we want to prove the convergence of the algorithm dSj) to minimizers of J. In 
order to do that, we need a characterization of solutions of the minimization problem ([3]) as 
the one provided in [381 Proposition 4.1] for the continuous setting. We specify the arguments 
in [38\ Proposition 4.1] for our discrete setting and we highlight the significant differences with 
respect to the continuous one. 
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5.1 Characterization of Solutions 

We make the following assumptions: 

{A^) (/5 : M — > M is a convex function, nondecreasing in M"*^ such that 

(i) (^(0) = 0. 

(ii) There exist c > and 6 > such that cz — b < ^p{z) < cz + b, for all z G M"*". 

The particular example we have in mind is simply (/^(s) = s, but we keep a more general 
notation for uniformity with respect to the continuous version in \38\ Proposition 4.1]. In this 
section we are concerned with the following more general minimization problem 

argmin„g^{:r^(n) := \\Tu - gg + 2a^(|Vtx|)(l^)} (18) 

where g £ 7i is a datum, a > is a fixed constant (in particular for (p{s) = s). 

To characterize the solution of the minimization problem ()18|) we use duality results from 
|22j . Therefore we recall the definition of the conjugate (or Legendre transform) of a function 
(for example see (22l Def. 4.1, pag. 17]): 

Definition 5.1. Let V and V* be two vector spaces placed in the duality by a bilinear pairing 
denoted by (•,•) and (j) : V ^ M be a convex function. The conjugate function (or Legendre 
transform) (f)* : V* ^ M. is defined by 

4>*{u*) = sup{{u,u*) — 4'{u)}. 

u€V 



Proposition 5.2. Let (^,u e Ti. If the assumption (A^) is fulfilled, then ( G djip{u) if and 
only if there exists M = {Mo,M) GfixTi'^, ^-^^ < ci G [0, +oo) for all x e Q such that 

{M{x),{Vu){x))^d +2a^{\{Vu){x)\) +2a(pl ( ^^^a^^ ) " ° for all x e n (19) 

T*Mo-dWM + C = (20) 
-Mo = 2{Tu-g), (21) 

where ipl is the conjugate function of (pi defined by ^i{s) = f{\s\), for s G M. 

If additionally ip is differentiable and |(Vn)(x)| ^ for x G f^, then we can compute M 

as 

The proof of this proposition specifies the one of |38l Proposition 4.1] to our discrete 
setting, it is technical, and it is deferred to the Appendix. 

Remark 5.3. (i) Forip{s) = s the function ipi from Proposition \5.2\ turns out to be <fi{s) = 
\s\. Its conjugate function ipl is then given by 



^p*i{s*) = sup{(s*,s) - 



for \s*\ < 1 
oo else 

Hence condition ()19p specifies as follows 

(M(x), (Vn)(x))K, + 2a\{Vu)ix)\ = 

and, directly from the proof of Proposition \ 5.S\ in the Appendix, \M{x)\ < 2a for all 

X en. 
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(ii) We want to highlight a few important differences with respect to the continuous case. 
Due to our definition of the gradient and its relationship with the divergence operator 
— div = V* no boundary conditions are needed. Therefore condition (10) of |22l Propo- 
sition 4-1] has no discrete correspondent in our setting. The continuous total variation 
of a function can be decomposed into an absolute continuous part with respect to the 
Lebesgue measure and a singular part, whereas no singular part appears in the discrete 
setting. Therefore condition (6) and (7) of Proposition 4-1] have not a discrete 
correspondent either. 

(Hi) An interesting consequence of Proposition [57B is that the map Sa = {I—PaK) is bounded, 
i.e., \\Saiz'')\\2 — ^ oo if and only if \\z''\\2 oo, for k — > oo. In fact, since 

Sa{z) = argmin — z||2 + 2a|Vu|(r2), 



from (j20p and (j2ip . we immediately obtain 



Sa{z) = z — - divM, 



and M is uniformly bounded. 
5.2 Convergence properties 

We return to the sequential algorithm ([5]). Let us explicitly express the algorithm as follows: 
Pick an initial Vi + V2 3 uf'^ + u^2^ := u^^^ € H, for example, uf^ = 0,i = 1,2, and iterate 

4"+i'^+i) = argmin (^i + nS"+^'')) £ = 0,...,L-1 

Tr|riUi=0 

(n+1,0) _ ~(n) 
^2 ~ ""2 

^(n+i,m+i) ^ a,rgmin ^^gv-^ J2iu^i^^^'^^ + U2, u^2^^^'"^^) m = 0, . . . , M - 1 
Tr|r2«2=0 
^(n+1) ^^^(n+l,L)^^{n+l,Af) 



~{n+l) 

~{n+l) 
k U2 



(23) 

Note that we do prescribe a finite number L and M of inner iterations for each subspace 
respectively and that n^""*"^-' = -u^"^"^^ + u'^~^^\ with u^"''^^^ 7^ u[^^^\ i = 1,2, in general. In 
this section we want to prove its convergence for any choice of L and M. 

Observe that, for a € and ||T|| < 1, 

\\ui — a\\2 — \\Tui — Ta\\2 > C\\ui — a|||, (24) 
for C = (1 - llTf) > 0. Hence 

J{u) = Ji{u,u^)<Ji{u,a), (25) 

and 

Jfiu, a) - Jf{u, u^) > C\\ui - a\\l (26) 
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Proposition 5.4 (Convergence properties). Let us assume that \\T\\ < 1. The algorithm in 
[23^) produces a sequence (u("^),igN '^^^h the following properties: 



(i) J(u(")) > J(u('^+i)) for all n en (unless u^") = n("+i) ); 
(ii) lim„^oo||n("+i) -u(")||2 = 0; 

(Hi) the sequence {u^'^'^)n&i has subsequences which converge in H. 
Proof. Let us first observe that 

By definition of u["^^'^^ and tlie minimal properties of u^^^'^^ in ()23p we liave 



From (1251) we liave 



Putting in line these inequalities we obtain 
In particular, from (j26p we have 
After L steps we conclude the estimate 
and 

j(^W) _ + 4")) > 114"+^''+'^ - 4"+'''^iii. 

By definition of ^2'^'''^'^^ and its minimal properties we have 

j(4"+^'^) +4")) > :ri(4"+^'^) +4"+^'^\4"+^'°)). 

By similar arguments as above we finally find the decreasing estimate 

> + 4"+^'"')) = J{u^^+'^) = Jiu^r^'^ + 4"+')), (27) 

and 

J(^/("V J(^x("+i)) 

\i=0 ■m=0 / 



OVERLAPPING DOMAIN DECOMPOSITION METHODS FOR TV-MINIMIZATION 15 



which verifies (i). 

From ( 1271 ) we have > By the coerciveness condition (C) is uni- 

formly bounded in Ti., hence there exists a convergent subsequence (n^"*^)^^^ and hence (iii) 
holds. Let us denote u^°°^ the limit of the subsequence. For simplicity, we rename such a sub- 
sequence by {u^"'^)n^n- Moreover, since the sequence {J'{u^"'^))n^n is monotonically decreasing 
and bounded from below by 0, it is also convergent. From (j28p and the latter convergence we 
deduce 



^L-l M-1 
^£=0 m=0 



+ E Il4"^''"^'^ - 4"^''"^ Hi ^ 0, n ^ oo. (29) 



In particular, by the standard inequality (a^ + 6^) > ^(a + ^)^ for a, 6 > and the triangle 
inequality, we have also 

II^W _^(n+i)||2 ^0, n^oo. (30) 
This gives (ii) and completes the proof. □ 
The use of the partition of unity {XI1X2} allows not only to guarantee the boundedness 

)nGN- 

Lemma 5.5. The sequences (£t['^'')„,gN o^n-d {u'^^)neN produced by the algorithm (|23p are 
bounded, i.e., there exists a constant C > such that Wu] II2 < C fori = 1,2. 



of (u*^"'^)„eN> but also of the sequences (n^"^)„gN and ({'■^"^'' 



Proof. From the boundedness of (u('^))„gN we have 



|uf)||2 = ||xin(")||2 < /^h(")||2 < C for i = 1,2. 



□ 



From Remark 15.31 (iii) we can also show the following auxiliary lemma. 

Lemma 5.6. The sequences (??i"'^^)n CLnd ^r]^'^^^)n o-re bounded. 
Proof. From previous considerations we know that 

4"'^-^) = 5«(z("'^'^"^) + 4"'^) - vt''^) - 4"'^). 

Assume {rj^'^^)n were unbounded, then by Remark [5. 31 (iii). also Sa{z^i^'^ +^2^ ~4"'^^) 
would be unbounded. Since {u'^^)n and (n^"''^'')„ are bounded by Lemma [5.51 and formula 
(p9]) . we have a contradiction. Thus (4"'^^)n has to be bounded. With the same argument 
we can show that (4"'*^^)" is bounded. □ 

We can eventually show the convergence of the algorithm to minimizers of J^. 

Theorem 5.7 (Convergence to minimizers). Assume \\T\\ < 1. Then accumulation points of 
the sequence (n("))„gp} produced by algorithm k23\) are minimizers of J . If J has a unique 
minimizer then the sequence (u^"^)„gN converges to it. 
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Proof. Let us denote u^°°^ the limit of a subsequence. For simplicity, we rename such a sub- 
sequence by (u("))„,gN. Prom Lemma 1^31 we know that ({t^"^)„gN> {u'^^)n£N and consequently 
{u^^'^^ )n£N ,{u^2^'^'^'')n£N a^e bounded. So the limit u^°°^ can be written as 

where is the limit of (u^"''^'*)„gN) ^2°°^ limit of (n2"'^^^)neN; and 'u^°°^ is the limit 

of for i = 1,2. Now we show that 7X2°°^ = By using the triangle inequality, 

from ([29j) it directly follows that 

||4"+i'*^)_4-)||2^0, n^oo. (32) 

Moreover, since X2 € is a fixed vector which is independent of n, we obtain from Proposition 
El (ii) that 

||^2(n(") -7x("+i))||2 ^0, n^oo, 

and hence 

112^0, n^oo. (33) 



Putting (I32p and (j33p together and noting that 

II (n+l,M) (n)|| , ||--{n) ^|i,(n+l,M) --(n+l),, 

11^2 ~ ^2 II2 "T 11^2 ~ ^2 II2 — 11^2 ~ ^2 II2 

we have 

||4-+M^)_4«+i) 112^0, n-00, (34) 

which means that the sequences (u2"'*^^)nGN and ({t2"^)ngN have the same limit, i.e., ^2°°^ = 
, which we denote by ^2°°^ . Then from and ([3T]) it directly follows that = u^^^ . 
As in the proof of the oblique thresholding theorem we set 



Fi(t.;"+^'^)) := ||ni"+^'^) - + 2a|V(ni"+^'^) + 4"^ )|(f^i) 

Six 

where 



The optimality condition for n^"^^'^^ is 



(n+l,L)x „ (n+l,L) 



0Gay,Fi(ni'™)+2??l 
where 

In order to use the characterization of elements in the subdifferential of |Vii|(il), i.e.. 
Proposition 15.21 we have to rewrite the minimization problem for Fi . More precisely, we 
define 



^^^^{n+i,L)^ := - 4"^ - 4"^''"'||2 + 2a|V(er^^'"0l(f^i) 



(n+l,i)||2 , o^,IV7^c("+l'^)^ 
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for g Vi with Tr jr^ = Then the optimahty condition for ig 

0G9A(ei"^''^^) + 27?i"+'''') (35) 
Note that indeed is optimal if and only if u^""*^^'^^ = ^j""'"^'^) _ is optimal. 



Analogously we define 



for g V2 with Tr = u^^~^^'^\ and the optimality condition for 

is 

G aACe^"""''*'^) + 24"+'-*') (36) 

where 

Let us recall that now we are considering functionals as in Proposition 15.21 with (p{s) = s, 
T = I, and Q = Qi, i = 1,2. Prom Proposition 15.21 and Remark 15.31 we get that ^ and 

consequently u^""*"^'^^ is optimal, i.e., —2r]^~^^'^^ € dFi{^^^^'^^), if and only if there exists 
an Mi^"+^^ = (M(5"+^\m{"+^^) e Vi x Vf with \m[''^^\x)\ < 2a for all x e fli such that 

(M^+^)(x),(V(^;"+^'^) +4")))(x))^, +2av,(|(V(n(-+^'^) +4")))(x^ = (37) 
-2(ni"+''^) (x) - 4"+'''') (x)) - div Af{"+') (x) - 27?^"+^'^) (x) = 0. (38) 

for all X G f^i. Analogously we get that ^2""*"^'^^ ^^'^ consequently u'^'^^'^^^ is optimal, i.e., 
_24"+i'*^) G 5F2(d"^'''''^)> if and only if there exists an iv4"+') = (M^;2+'\ M^"+')) G 
^2 X ^2'' with |M^"+^^(x)| < 2q for all x e ^2 such that 

(M("-^^)(x),(V(nr''^)+4"+^'^'^)))(x))^.+2a^(|(V(4"+^'"^ (39) 
-2(4"+^'*^) (x) - 4"+''"')(x)) - divMf +^)(x) - 24"+^'^)(x) = 0, (40) 

for all X G 0,2- Since {M["'\x))neN is bounded for all x G fii and (M2"^(^))nGN is bounded 
for all X G ^2, there exist convergent subsequences {M["'^\x))k(^N and (M2"'°^(^))fceN- Let 
us denote M-[°°^(x) and M2°°''(x) the respective limits of the sequences. For simplicity we 
rename such sequences by {M["\x))neN and (M2"^(^))neN- 

Note that, by Lemma [5^61 (or simply from ([38|) and (00])) the sequences (77i"''^^)nGN and 
(4"'^^^)neN are also bounded. Hence there exist convergent subsequences which we denote, 
for simplicity, again by (f?["''^^)neN and (r/2"'^'*)nGN with limits r]^'^\ i = 1,2. By taking in 
([571) - pn|) the limits for n — 00 we obtain 

{m[^\x), (V(nS~) + 4°"^))(x))iRd + 2a(^(|(V(nS~) + 4"^))(x)|) = for ah x G f^i 
-2{u'f°^ (x) - (x)) - div m{°°) (x) - 24°°) (x) = for ah x G f^i 



>+l)||2 , o^,|v7/t(n+l,A/)^ 



OVERLAPPING DOMAIN DECOMPOSITION METHODS FOR TV-MINIMIZATION 18 



{M^°°\x), (V(ni~^ + + 2av9(|(V(ni~^ + = for all x € O2 

-2(4°°) (x) - z^~) (x)) - div M2^~) (x) - 2^^) (x) = for all x £ ^2 

Since suppry^""-* = Ti and supprj^^ = T2 we have 

{m[°°\x), (V(n(°°))(x))Kd + 2a(^(|(V(n(°°))(x)|) = for all x G J7i 

-2r*((ru(°°))(x) - 5(°°)(x)) - divM{°°)(x) = for all x G i7i \ Ti 

(m(°°)(x), (V(n(~))(x))i,. + 2a(^(|(V(n(~))(x)|) = for all x G ^^^^ 
-2r*((rn(~))(x) - c/(°°)(x)) -divM^°°)(x) = for all x G Jl2\r2. 

Observe now that from Proposition 15.21 we also have that G J'{u^°°^) if and only if there 
exists M(°°) = (M^°°\m(°°)) with \M^°°\x)\ < 2a for all x G such that 

(M(°°)(x), (V('u(°°))(x))Kd + 2a(/j(j(V('u(°°))(x)|) = for ah x G !^ 

-2r*((ru(°°))(x) - ^(^^(x)) - divM(°°)(x) = for ah x G ^^^^ 



Note that M^°°\x),j = 1,2, for x G r^i n O2 satisfies both (jH]) and (jH]). Hence let us 
choose 

[m{°°)(x) ifxGf^i\ri 
]m^°°)(x) if X G (1^2 UTi 



X 



With this choice of M^°°^ equations ([^T|) - (jl3]) are valid and hence n(°°) is optimal in ri. □ 

Remark 5.8. (i) If'Vu^°°\x) 7^ for x G Oj, j = 1, 2, then Mj°°^ is given as in equation 
(EID by 

M(~)(x)--2ai^^^i!lM^ 
^ '"|(V^Mb,)(x)r 

(^iij T/ie boundedness of the sequences (u^^)n£N o-nd (u2'^))„gN has been technically used for 
showing the existence of an optimal decomposition + ^2°°) in the proof 

of Theorem \5. 7\ Their boundedness is guaranteed as in Lemma 15.51 by the use of the 
partition of the unity {xitX2}- Let us emphasize that there is no way of obtaining the 
boundedness of the local sequences {u^-P'^^)n£N o-nd {u'^'^^^)n£N otherwise. In Figure\^ 
we show that the local sequences can become unbounded in case we do not modify them 
by means of the partition of the unity. 

( Hi ) Note that for deriving the optimality condition (|43p for u^°°^ we combined the respective 
conditions ()4ip and (I42p for and u'^^ . In doing that, we strongly took advantage 
of the overlapping property of the subdomains, hence avoiding a fine analysis of r/|°°^ 
and T]^^ on the interfaces Fi and T2 . This is the major advantage of this analysis with 



respect to the one provided in 126^ for nonoverlapping domain decompositions. 
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6 A parallel algorithm and its convergence 

The parallel version of the previous algorithm (j23p reads as follows: Pick an initial V1 + V2 B 



~(o) , ~(o) 



€ TC, for example uf^^ = 0, i = 1, 2, and iterate 



u 



u 



u 



(n+1,0) _ -in) 
1 — "' 

(n+ll+l) 
1 



argmin „^gvi J{{ui + u''^\u^^~^^'''') £ = 0,...,L-1 



(n+1,0) ~{n) 
U2 = U2 

(n+l,m+l) 



Trlri«i=0 



Uo ' = argmm J2[Ui +U2.,u\ ') m = 0, . . . , M — 1 

Tr|r2M2=0 



(n+l,m) 



(44) 



(n+l) 
r,{"+l) 



X2 • 



We are going to propose similar convergence results as for the sequential algorithm. 

Proposition 6.1 (Convergence properties). Let us assume that \\T\\ < 1. The parallel algo- 
rithm ^4\ ) produces a sequence (ii^"'')„,gN T~(- with the following properties: 

(i) J'(u(")) > J^(u("+i)) for all n en (unless u^") = u("+^) ); 

(ii) lim„^oo||n("+^) -nW||2 =0; 

(Hi) the sequence (ti''"'')nGN has subsequences which converge in Ti. 
Proof. With the same argument as in the proof of Theorem 15.41 we obtain 

L-l 



j(^(-)) _ J(^V' 



(n+l,L) ~(n) 



+ u\ 



> 



(n+l,e+l) (n+l/) 1,2 



1=0 



and 



J{n(^))-J{u^;^> +u)2 



(n) (n+l,M)s 



M-1 



> 



(n+l,m+l) (n+l,m)||2 



tin 



m=0 



Hence, by summing and halving 



> 



C 



1 



(n+l,L) ~(n) 



+ u^"0 + ^(^'r+«2 



(n) (n+l,M) 



)) 



(n+l,£+l) (n+l/)||2 



+ E 



I (n+l,m+l) (n+l,m)|i2 
— ^2 II2 



m=0 



We recah that J(u('^)) = ||rti(") -g\\l + 2a\Vu^'^^\{n) (and T is linear). Then, by the standard 
inequality (a^ -\-b'^) > |(a + b)"^ for a, 6 > 0, we have 



T 



(n+l,L) , (n+l,M) 



) + ^ 



< i||T(tzi"+''^)+4")) 



.lli + ^l|T(4"^ + nr^'*^V^^Ili- 
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Moreover we have 

By the last two inequahties we immediately show that 
hence 



C 
> — 
- 2 



^L-l M-1 \ 

^||^(n+l,m)_^(n+l,.)||2+ ^ ||^(n+l,m+l)_^(„,+l,m)||2 (45) 



v^=0 m=0 



Since the sequence {J'{u^'^^))neN is monotonically decreasing and bounded from below by 0, 
it is also convergent. From (|45|) and the latter convergence we deduce 



/L-l M-i \ 

E IKS"^'''^'^ - 4^^^'''^ IIh + E Il4"^'''"^'^ - 4"^'''"^ Hi - 0, n - oo. (46) 

\£=0 m=0 / 

In particular, by again using (a^ + 6^) > ^(fl + for a, 6 > and the triangle inequality, we 
also have 

II^W _^('^+i)||2 ^ 0, n^oo. (47) 
The rest of the proof follows analogous arguments as in that of Proposition 15. 4[ □ 

Analogous results as the one stated in Lemma [5. 51 and Lemma [5. 6 1 also hold in the parallel 
case. With these preliminary results the following theorem holds: 

Theorem 6.2 (Convergence to minimizers). Assume \\T\\ < 1. Then accumulation points of 
the sequence (n("))„gp} produced by algorithm ^44\ ) minimizers of J . If J has a unique 
minimizer then the sequence (u^"^)nGN converges to it. 

Proof. Note that u^""*"^-* is the average of the current iteration and the previous, i.e., 

(n+l,L) {n+l,M) , „,(n) 



Observe that the sequences (u[^"^^'^^)„gN5 (4"~'"^'^^^)"GN and (n'^"))„gN are bounded. Hence 
there exist convergent subsequences. By taking the limit for n ^ cxd we obtain 



which is equivalent to 

(oo) (oo) , (oo) 



With this observation the rest of the proof follows analogous arguments as in that of Theorem 
EH □ 
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7 Applications and Numerics 

In this section we shall present the application of the sequential algorithm ([5]) for the min- 
imization of J in one and two dimensions. In particular, we show how to implement the 
dual method of Chambolle [10] in order to compute the orthogonal projection PaK{g) in the 
oblique thresholding, and we give a detailed explanation of the domain decompositions used 
in the numerics. Furthermore we present numerical examples for image inpainting, i.e., the 
recovery of missing parts of images by minimal total variation interpolation, and compressed 
sensing [H [71 El |21], the nonadaptive compressed acquisition of images for a classical toy 
problem inspired by magnetic resonance imaging (MRI) [71 [29] . The numerical examples of 
this section and respective Matlab codes can be found at |40j . 

7.1 Computation of PaK{g) 

To solve the subiterations in ([5|) we compute the minimizer by means of oblique thresholding. 
More precisely, let us denote U2 = u'^\ ui = u^^~^^'^~^^\ and zi = u^^'^^'^'^ +T*{g — Tu2 — 
Tu^~^^'^^). We shall compute the minimizer ui of the first subminimization problem by 

«! = (/- PaK){zi + U2 - V) - U2 ^Vi 

for an r/ G Vi with suppr/ = Ti which fulfills 

Tr (r/) = Tr |ri (zi + Paxiv - zi - U2)) . 
Hence the element € Vi is a limit of the corresponding fixed point iteration 

r/(0) G yi,suppr?W = n, r?(™+i) = (Tr |rj* Tr |r, (^i + PaKiv^""^ - zi - U2)) , m > 0. 

(48) 

Here K is defined as in Section [21 i.e., 

K= |divp:pG W*^, |p(x)|^ < 1 VxGr?}. 

To compute the projection onto aK in the oblique thresholding we use an algorithm proposed 
by Chambolle in [10]. His algorithm is based on considerations of the convex conjugate of 
the total variation and on exploiting the corresponding optimality condition. It amounts to 
compute Paxig) approximately by adivp*^""^, where p*-""^ is the nth iterate of the following 
semi-implicit gradient descent algorithm: 

Choose r > 0, let p^^^ = and, for any n > 0, iterate 

in+i)(. ^ p(")(rr)+r(V(divp(")-ff/a))(:r) 
^' 1 + r |(V(divp(") -5/a))(x)| 

For r > sufficiently small, i.e., r < 1/8, the iteration adivp^"^ was shown to converge 
to PaK{g) as n —f 00 (compare [lOl Theorem 3.1]). Let us stress that we propose here this 
algorithm just for the ease of its presentation; its choice for the approximation of projections 
is of course by no means a restriction and one may want to implement other recent, and 
perhaps faster strategies, e.g., [TT l [TUl [271 1^ 
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7.2 Domain decompositions 

In one dimension the domain = [a, b] is split into two overlapping intervals. Let nO,2\ =: 
G be the size of the overlap of Oi and $72- Then we set =: ni = \ ^2^ ~\ , = 
and ^2 = [ni — G + 1, b]. The interfaces Fi and r2 are located in i = ni + 1 and ni — G 
respectively (cf. Figure [2]). The auxiliary functions xi X2 can be chosen in the following 
way (cf. Figured]): 

1 Xj G \ ^2 

1 - ^(i - (ni - G + 1)) Xieninn2 

1 Xj G r22 \ ^1 

^(i - (ni - G + 1)) Xieninn2' 

Note that Xi(^i) + X2{xi) = 1 for all Xj € (i.e for all i = 1, . . . , N). 



Xiixi) = 
X2{Xi) = 



2 




Figure 1: Auxiliary functions xi a-nd X2 for a-n overlapping domain decomposition with two subdo- 
mains. 

In two dimensions the domain Vl = [a, b] x [c, d\ is split in an analogous way with respect 
to its rows. In particular we have = [a, ni] x [c, d\ and O2 = [ni — G + 1, b\ x [c, d], compare 
Figure [31 The splitting in more than two domains is done similarly: 

Set = r^iU. . .ur2_A/, the domain $7 decomposed into M domains Oj, z = 1, . . . ,7V, 
where i7j and fij+i are overlapping for i = l,...,7V — 1. Let fi =: G 

equidistant for every i = l,...,7V— 1. Set s = \Ni/J\r\ . Then 

J7i = + X [c,d] 
for i = 2 : TV - 1 

G G 

!7i = [{i - l)s - - + 1, is + -] X [c, d] 

end 

= [(A/--l)s-^ + l,Ar^] X [c,d]. 
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The auxiliary functions Xi can be chosen in an analogous way as in the one dimensional case: 
'^(ii-((i-l)s-G/2 + l)) ixi,,yi,) eni^iDQi 



for i = 1, . . . , A/" with JIq = ^A^+i 



l-^{h-{is-G/2 + l)) {xi^,yi2) G^nf^i+i 



r2 

-e — 



Hi 

Figure 2: Overlapping domain decomposition in ID. 



a = xi 


ni\n2 


Xni-G 











Ti 


b = Xjy 





Figure 3: Decomposition of the image in two domains Qi and 



To compute the fixed point rj of (jlSp in an efficient way we make the following considera- 
tions, which allow to restrict the computation from to a relatively small stripe around the 
interface. The fixed point r] is actually supported on Fi only, i.e., rj(x) = in f^i \Fi. Hence, 
we restrict the fixed point iteration for to a relatively small stripe f^i C l^i Analogously, one 
implements the minimizations of 772 on 0,2- A similar trick was also used in [26] to compute 
suitable Lagrange multipliers at the interfaces of the nonoverlapping domains. However, there 
we needed to consider larger "bilateral stripes" around the support of the multiplier, making 
the numerical computation slightly more demanding for that algorithm. 

7.3 Numerical experiments 

In the following we present numerical examples for the sequential algorithm (|23p in two 
particular applications: signal interpolation/image inpainting, and compressed sensing. 

In Figure H] and Figure [5] we show a partially corrupted ID signal on an interval 0, of 
100 sampling points, with a loss of information on an interval D d Vl. The domain D 
of the missing signal points is marked with green. These signal points are reconstructed 
by total variation interpolation, i.e., minimizing the functional ^7 in ([3]) with a = 0.4 and 
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overlap 5 






" overlap 1 























Figure 4: We present a numerical experiment related to the interpolation of a ID signal by total 
variation minimization. The original signal is only provided outside of the green subinterval. The 
initial datum g is shown in (a). As expected, the minimizer is the constant vector 1, as shown in 
(b). In (c) and (d) we display the rates of decay of the relative error and of the value of J respectively, 
for applications of the algorithm with different sizes G=l,5,10,20,30 of the overlapping region of 
two subintervals. 
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Figure 5: We show a second example of total variation interpolation in ID. The initial datum g is 
shown in (a). As expected, a minimizer is (nearly) a piecewise linear function, as shown in (b). 
In (c) and (d) we display the rates of decay of the relative error and of the value of J' respectively, for 
applications of the algorithm with different sizes G=l,5,10,20,30 of the overlapping region of two 
subintervals. 
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Tu 



1 



n\D 



where Iq\d is the indicator function of 0,\ D. A minimizer u^°°^ of J' is 
precomputed with an algorithm working on the whole interval 0, without any decomposition. 
We show also the decay of relative error and of the value of the energy J' for applications of 
algorithm (j23p on two subdomains and with different overlap sizes G = 1,5, 10,20,30. The 
fixed points t^'s are computed on a small interval Oj, i = 1,2, of size 2. These results confirm 
the behavior of the algorithm (j23p as predicted by the theory; the algorithm monotonically 
decreases and computes a minimizer, independently of the size of the overlapping region. 
A larger overlapping region does not necessarily imply a slower convergence. In these figures 
we do compare the speed in terms of CPU time. In Figure [6] we also illustrate the effect 
of implementing the BUPU within the domain decomposition algorithm. In this case, with 
datum g as in Figure [Sj we chose a = 1 and an overlap of size G = 10. The fixed points r/'s 
are computed on a small interval Oj, i = 1,2 respectively, of size 6. 



1 00 iterations 



Reconstruction 
O Interface 



- u2 

Inpainting Region 





Reconstruction 
O Interface 



- u2 

Inpainting Region 



(a) 



(b) 



Figure 6: Here we present two numerical experiments related to the interpolation of a ID signal by 
total variation minimization. The original signal is only provided outside of the green subinterval. On 
the left we show an application of algorithm (j23p when no correction with the partition of unity is 
provided. In this case, the sequence of the local iterations u^^\u^'^ is unbounded. On the right we 
show an application of algorithm (|23p with the use of the partition of unity which enforces the uniform 
boundedness of the local iterations u^^\ u'^'^ . 

Figure [7] shows an example of the domain decomposition algorithm (j23p for total variation 
inpainting. As for the ID example in Figures [4]l6] the operator T is a multiplier, i.e., Tu = 
^n\D ' "^^j where ft denotes the rectangular image domain and D C ^ the missing domain in 
which the original image content got lost. The regularization parameter a is fixed at the value 
10~^. In Figure [7] the missing domain D is the black writing which covers parts of the image. 
Here, the image domain of size 449 x 570 pixels is split into five overlapping subdomains with 
an overlap size G = 28 x 570. Further, the fixed points r^'s are computed on a small stripe 
r^j, i = 1, . . . , 5 respectively, of size 6 x 570 pixels. 

Finally, in Figure [8] we illustrate the successful application of our domain decomposition 
algorithm (j23p for a compressed sensing problem. Here, we consider a medical-type image 
(the so-called Logan-Shepp phantom) and its reconstruction from only partial Fourier data. 
In this case the linear operator T = S oj^, where J-" denotes the 2D Fourier matrix and 5 is a 
downsampling operator which selects only a few frequencies as output. We minimize J' with 
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(a) 



(b) 



Figure 7: This figure shows an application of algorithm (|23|) for image inpainting. In this simulation 
the problem was split into five subproblems on overlapping subdomains. 



Q set at 0.4 x 10~^. In the application of algorithm (j23p the image domain of size 256 x 256 
pixels is split into four overlapping subdomains with an overlap size G = 20 x 256. The fixed 
points r/'s are computed in a small stripe Oj, i = 1, . . . , 4 respectively, of size 6 x 256 pixels. 
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A Proof of Proposition 15.21 

It is clear that C ^ dj'^{u) if and only if u = argmin„g^{j7'^(?;) 
the following variational problem: 



ml{J^{v) - {C,v)n} = inf {||Tz; -^^11^ + 2a^{\Vvm) - {C,v)n} 



(C) ^)'h}j and let us consider 

{V) 
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Sampling domain in the frequency plane 





(a) 



(b) 





(c) 



(d) 



Figure 8: We show an application of algorithm ((23l) in a classical compressed sensing problem for 
recovering piecewise constant medical-type images from given partial Fourier data. In this simulation 
the problem was split via decomposition into four overlapping subdomains. On the top-left figure, 
we show the sampling data of the image in the Fourier domain. On the top-right the back-projection 
provided by the sampled frequency data together with the highlighted partition of the physical domain 
into four subdomains is shown. The bottom figures present intermediate iterations of the algorithm, 
i.e., ^(26) and u^^^s)^ 



OVERLAPPING DOMAIN DECOMPOSITION METHODS FOR TV-MINIMIZATION 29 



We denote such an infimum by inf • Now we compute (\V*\i the dual of . Let J- : Ti 
g :n xW^ ^R, Gi -.n ^R, 02 ■■H'^ such that 

J^{v) = -{C,v)h 
Gi{wo) = ||'"^o-5ll2 
g2{tB) = 2aip{\tjd\){n) 

g{w) = giiwo)-\-g2{w) 

with w = {wo,w) € 7^ X H'^. Then the dual problem of is given by (cf. [22l p 60]) 

sup {-T*{A*p*)-g*{-p*)} {V* 

where A : H ^ TC x is defined by 

Av = {Tv,{Vv)^,...,{Vv)'^) 



and A* is its adjoint. We denote the supremum in ([23) by sup (|'P*p . Using the definition of 
the conjugate function we compute JT* and g* . In particular 

r'{A*p*) = snp{{A*p*,v)n - Hv)} = sup(AV + C,v)n = C ^5 " " 

veH v&H loo otherwise 

where p* = (pqjP*) ^'^d 

g*{p*)= sup {{j)*,w)j^y,yid-g{w)} 

= sup {{vl,WQ)ri + {p*,w)^d-gi{wQ)-g2{w)} 

w=[wo,u))e'Hx'H'^ 

= sup {{pl,wo)n - gi{wo)} + snp {{p*,w)f^d - g2{w)} 



We have that 



gi{p*) = 2a^l['^^]m 



and (see [22] 

G*(n*\ = Imn* ( 

^ 2a 

if ^ Domipl, where ip'^ is the conjugate function of (pi defined by 

^i(s) :=ip{\s\) s€M. 

For ease we include in Appendix [B| the explicit computation of these conjugate functions. So 
we can write in the following way 
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where 

JC = !^p* e n X : G Bom^l for all x e n, A*p* + ( = o| . 

The function fi also fulfills assumption (^^)(ii) (i.e., there exists ci > 0, 6 > such that 
ciz — b < ^i{z) < ciz + b, for all z € M"*"). The conjugate function of ipi is given by 
(/9*(s) = sup^g]g{(s, z) — ifi{z)}. Using the previous inequalities and that ipi is even (i.e., 
(/9i(z) = ipi(—z) for all z € M) we have 



(sup{(s, z) — ci\z\ + 6} >) sup{(s, z) — ^fiiz)} > sup{(s, z) — ci\z\ — b} = 
zeM. 2gM zeM. loo else 

In particular, one can see that s € Domipl if and only if |s| < ci. 
From A*p* + = we obtain 



b if |s| < ci 
(1.50) 



{A*p*,uj)'H + {C,uj)n = {p*,Auj)y^d+i + {C,uj)'H = {po,Tuj)n + {p*,Vuj)-^d + (C,w>w = for all w G 
Then, since {p* ,'Vuj)i-^d = {— div p* ,uj)-}i (see Section[2]), we have 

T*p*Q -divp* + C = 0. 
Hence we can write /C in the following way 

]C= \ p* = {p*n,p*) eUxH'' : < d for all x e n, T*pl - divp* + C = 

[ 2a 

We now apply the duality results from [221 Theorem III. 4.1], since the functional in (fPl) is 
convex, continuous with respect to Av in Ti.xTi'^, and inf (jP]) is finite. Then inf (fP|) = sup (|'P*P € 
M and has a solution M = {Mq, M) E /C. 

Let us assume that u is a solution of (j^ and M is a solution of From inf (j^ = 

sup (|P*|) we get 



|Tn-5||^ + 2av?(|Vn|)(17)-(C,n)>^ = -(-^ + g,-Mo\ - 2a(^]: f ^ ) (17) (1.51) 



H 



2a 



where M = (Mo,M) G H x W^, < ci and r*Mo - divM + C = 0, which verifies the 

direct implication of (|20p . In particular 

-(C,n)H = {T*Mo,u)n - {dwM,u)n = {Mo,Tu)n + (M, Vn)^^, 

and 

\\Tu-g\\l+{Mo,Tu)n+{M,Vu)^d+2a^{\Vum)+(^^ + g,-Mo'^ +2a^l (^^^ (0) = 0. 

(1.52) 

Let us write (I1.52P again in the following form 

d 

\{Tu - g){x)\' + Mo{x){Tu){x) + ^ E ^'■(^)(Vu)^'(^) + E 2«¥'(l(Vtx)(x)|) 
x€n xeQ j=i xen 

+ E (^ + .w) (-MoW) + E (%^) = »■ 

(1.53) 

Now we have 
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1. 2a^(|(Vn)(x)|) + E?=iM^'(^)(Vny(x) + 2a(^t (MM) > 2a^(|(Vn)(x)|) 

E?=i |M^(x)||(Vn)-' (x)|+2a^^ (^^) > by the definition of (^^ , since 2099^ 
sup5gK4(iVP(x),5)i,. -2a(^(|5|)} =sup5gjj4(|M^-(x)|,|5|)K. -2a(^(|5|)} . 

2. KT^x - g){x)\^ + Mo(x)(r^.)(x) + + 5(x))(-Afo(x)) = (((Tn)(x) - ^^(x)))^ + 
Mo{x){{Tu)ix) - g{x)) + (^)' = [{iTu)ix) - g{x)) + > 0. 

Hence condition (jl.52p reduces to 

2a(/p(|(Vn)(x)|) +^M^(x)(Vn)^(x) + 2a^t 0-^^) = for all x € !^ (1.54) 
j=i 

- Mo(x) = 2{{Tu){x) - g{x)) for all x G 0. (1.55) 

Conversely, if such an M = (Mo,M) € Ti x with < ci exists which fulfills 

conditions p^ - (j2T]) . it is clear from previous considerations that equation (|1.5ip holds. Let 
us denote the functional on the left side of (ll.51|) by 

P(n) := \\Tu - g\\l + 2a^(|Vtx|)(0) - (C, u)h 

and the functional on the right side of (11.511) by 

P*{M) ■.= -(zMl + g^ -Mo^^ - 2a^\ (J7). 



We know that the functional P is the functional of (|^ and P* is the functional of (|'P*p . 
Hence inf P = inf (lP]l and supP* = sup (IP* p . Since P is convex, continuous with respect to 
kum.'H-K Ti'^, and infdEj) is finite we know from duality results |22l Theorem III. 4.1] that 
inf(l£l)= sup(l2!])e M. We assume that M is no solution of i.e., P*{M) < supdH]), and 

« is no solution of ([P|) . i.e, P(u) > inf (jPj) . Then we have that 



P{u) > inf dB = sup ^ > P*{M). 

Thus (jl.Sip is valid if and only if M is a solution of (|'P*P and it is a solution of ([2]) which 
amounts to saying that ( G dj^p{u). 

If additionally is differentiable and |(Vn)(x)| 7^ for a; G $7, we show that we can 
compute M{x) explicitly. From equation (I19p (resp. (jl.54p ) we have 

2a^\ (^1^) = -(^(^)' (V^)(^))m^ - 2a^(|(Vn)(x)|). (1.56) 
From the definition of conjugate function we have 



2asup<! ( -(^i(t) 



t>o 



2asupsup|( — ,3) -ipi{\S\) 

t>0 seRd I \ / 
\S\=t 

sup {{-M{x),S)^,-2a^{\Sm)}. 



(1.57) 
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Now, if |(Vm)(x)| 7^ for x G 0, then it follows from (jl.56p that the supremum is taken on 
in 5 = |(Vu)(x)| and we have 

Vsi-{M{x),S)^, - 2a^i\Sm)) = 

which implies 

|(Vn)(x)| 

and verifies (j22|) . This finishes the proof. 

B Computation of conjugate functions 

Let us calculate the conjugate function of the convex function Qi{wo) = \\wo — g\^. From 
Definition 15 . 1 1 we have 

Q\{vl)= sup {(u;o,pS)w - ^ilit'o)} = sup {(wo,^'o)w - (^0 - 5,u^o - S-)?^}- 

We set H{wo) ■= {woiPajn — {wo — g,WQ — g)'n. To get the maximum of H we calculate the 
Gateaux-differential at wq of H, 

H'iwo) =p*o- 2{wo -g) = 

and we set it to zero H'{wq) = 0, since H"{wo) < 0, and we get wq = ^ + g. Thus we have 
that 

sup H{wo) = (4 + 9,Po) =QI{PI) 



Now we are going to calculate the conjugate function of G2i'w) = 2aip{\w\){Q). Associated 

+ — TI3' + 



to our notations we define the space -H+ = R+^^''-''^^ From Definitional] we have 



02{P*) = sup {{w,p*)nd - 2af{\w\){n)} 

= sup sup {{w,p*)'y^d — 2aip{\w\){^)} 

\w{x)\=t{x) 

= sup{{t,\p*\)n-2aip{t){n)}. 



If ip were an even function then 



sup {{t, \p*\)n - 2a^{t)in)} = sup{(i, \p*\)n - 2a^{t){n)} 

= 2asup|/i,|^\ -^mn) 
ten l\ 2a 



where (p* is the conjugate function of 
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Unfortunately (f is not even in general. To overcome this difficulty we have to choose a 
function which is equal to ip{s) for s > and does not change the supremum for s < 0. For 
instance, one can choose ^pi{s) = ^{\s\) for s € M. Then we have 



where is the conjugate function of ipi. Note that one can also choose ^pi{s) = ^p{s) for 
s > and (pi{s) = 00 for s < 0. 
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